rm(list=ls())
library(DirectEffects)
library(readstata13)
library(lfe)

setwd("C:/Users/m1049/Dropbox/CPS Replication Final")

# road
df<-read.dta13("acde.dta")

m1<-felm(road ~ treatment_binary | id + year | 0 | id, 
         data = df)
summary(m1)

# bootstrap

set.seed(12345)
reps=1000
par.est<-matrix(NA,nrow=reps,ncol=1)
colnames(par.est)<-c("acde_pop_mean")


for (i in 1:reps){
  
  temp <- df[sample(1:nrow(df), replace=TRUE),]
  
  form_main <- road ~ treatment_binary + as.factor(id) +as.factor(year) | I(log(transfer)) | I(log(pop))+mean
  
  direct <- sequential_g(form_main, data = temp)
  
  par.est[i,1]<-direct$coefficients[2]
  
}

# baseline 6.148 

boot<-as.data.frame(par.est)

mean(boot$acde_pop_mean) # 6.24271

se<-sd(boot$acde_pop_mean); se

t=mean(boot$acde_pop_mean)/se; t

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.96*rep(se,each=2) 


##  water

m1<-felm(watersupply ~ treatment_binary | id + year | 0 | id, 
         data = df)

set.seed(12345)
reps=1000
par.est<-matrix(NA,nrow=reps,ncol=1)
colnames(par.est)<-c("acde_pop_mean")


for (i in 1:reps){
  
  temp <- df[sample(1:nrow(df), replace=TRUE),]
  
  # population + mean household income
  
  form_main <- watersupply ~ treatment_binary + as.factor(id) +as.factor(year) | I(log(transfer)) | I(log(pop))+mean
  direct <- sequential_g(form_main, data = temp)
  par.est[i,1]<-direct$coefficients[2]
  
}

boot<-as.data.frame(par.est)

mean(boot$acde_pop_mean) 

se<-sd(boot$acde_pop_mean); se

t=mean(boot$acde_pop_mean)/se; t

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.96*rep(se,each=2)

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.645*rep(se,each=2) 

# community

m1<-felm(lib ~ treatment_binary | id + year | 0 | id, 
         data = df)
summary(m1)

set.seed(12345)
reps=1000
par.est<-matrix(NA,nrow=reps,ncol=1)
colnames(par.est)<-c("acde_pop_mean")

for (i in 1:reps){
  
  temp <- df[sample(1:nrow(df), replace=TRUE),]
  
  # population + mean household income
  
  form_main <- lib ~ treatment_binary + as.factor(id) +as.factor(year) | I(log(transfer)) | I(log(pop))+mean
  direct <- sequential_g(form_main, data = temp)
  par.est[i,1]<-direct$coefficients[2]
  
}

# baseline -.5240538

boot<-as.data.frame(par.est)

mean(boot$acde_pop_mean) 

se<-sd(boot$acde_pop_mean); se

t=mean(boot$acde_pop_mean)/se; t

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.96*rep(se,each=2) 

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.645*rep(se,each=2) 

# center

m1<-felm(center ~ treatment_binary | id + year | 0 | id, 
         data = df)

summary(m1)

set.seed(12345)
reps=1000
par.est<-matrix(NA,nrow=reps,ncol=1)
colnames(par.est)<-c("acde_pop_mean")

system.time(for (i in 1:reps){
  
  temp <- df[sample(1:nrow(df), replace=TRUE),]
  
  # population + mean household income
  
  form_main <- center ~ treatment_binary + as.factor(id) +as.factor(year) | I(log(transfer)) | I(log(pop))+mean
  direct <- sequential_g(form_main, data = temp)
  par.est[i,1]<-direct$coefficients[2]
  
}
)

# baseline -1.166163  

boot<-as.data.frame(par.est)

mean(boot$acde_pop_mean) 

se<-sd(boot$acde_pop_mean); se

t=mean(boot$acde_pop_mean)/se; t

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.96*rep(se,each=2) 

# storage

m5<-felm(max ~ treatment_binary | id + year | 0 | id, 
         data = df)

summary(m5)

set.seed(12345)
reps=1000
par.est<-matrix(NA,nrow=reps,ncol=1)
colnames(par.est)<-c("acde_pop_mean")

for (i in 1:reps){
  
  temp <- df[sample(1:nrow(df), replace=TRUE),]
  
  # population + mean household income
  
  form_main <- max ~ treatment_binary + as.factor(id) +as.factor(year) | I(log(transfer)) | I(log(pop))+mean
  direct <- sequential_g(form_main, data = temp)
  par.est[i,1]<-direct$coefficients[2]

}

# baseline -184

boot<-as.data.frame(par.est)

mean(boot$acde_pop_mean) 

se<-sd(boot$acde_pop_mean); se

t=mean(boot$acde_pop_mean)/se; t 

rep(mean(boot$acde_pop_mean),each=2)+c(-1,1)*1.96*rep(se,each=2) 
